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NOTES  AND  CORRESPONDENCE 

Trajectory  Calculations  for  Spherical  Geodesic  Grids  in  Cartesian  Space 

Francis  X.  Giraldo 

Naval  Research  Laboratory,  Monterey,  California 


9  April  1998  and  30  July  1998 
ABSTRACT 

This  paper  shows  how  to  obtain  accurate  and  efficient  trajectory  calculations  for  spherical  geodesic  grids  in 
Cartesian  space.  Determination  of  the  departure  points  is  essential  to  characteristic-based  methods  that  trace  the 
value  of  a  function  to  the  foot  of  the  characteristics  and  then  either  integrate  or  interpolate  at  this  location.  In 
this  paper,  the  departure  points  are  all  computed  in  relation  to  the  spherical  geodesic  grids  that  are  composed 
of  a  disjoint  set  of  unstructured  equilateral  triangles.  Interpolating  and  noninterpolating  trajectory  calculation 
approaches  are  both  illustrated  and  the  accuracy  of  both  methods  are  compared.  The  noninterpolating  method 
of  McGregor  results  in  the  most  accurate  trajectories.  The  challenge  in  using  McGregor’s  method  on  unstructured 
triangular  grids  lies  in  the  computation  of  the  derivatives  required  in  the  high-order  terms  of  the  Taylor  series 
expansion.  This  paper  extends  McGregor’s  method  to  unstructured  triangular  grids  by  describing  an  accurate 
and  efficient  method  for  constructing  the  derivatives  in  an  element  by  element  approach  typical  of  finite  element 
methods.  An  order  of  accuracy  analysis  reveals  that  these  numerical  derivatives  are  second-order  accurate. 


1.  Introduction 

The  solution  of  partial  differential  equations  on  the 
sphere  is  of  prime  importance  in  meteorology  and 
oceanography.  The  proper  coordinate  system  would  ap¬ 
pear  to  be  the  spherical  coordinates  but  this  system  poses 
some  challenging  problems  at  the  poles  not  only  for 
Eulerian  formulations  of  the  equations  but  for  Lagrang- 
ian  formulations  as  well.  The  pole  problem  can  be  over¬ 
come  in  a  variety  of  ways,  such  as  the  use  of  Cartesian 
rather  than  spherical  coordinates  to  write  the  differential 
equations.  This  approach  was  used  in  Giraldo  (1997) 
for  the  advection  equation.  This  approach  may  also  be 
used  for  the  shallow  water  equations  by  using  the  La¬ 
grange  multiplier  formulation  of  Cote  (1988).  Cote 
writes  the  equations  in  Cartesian  coordinates  but  then 
includes  an  extra  forcing  term  obtained  by  using  La¬ 
grange  multipliers.  This  forcing  term  constrains  the  mo¬ 
tion  of  all  fluid  particles  to  remain  on  the  sphere.  We 
are  very  much  interested  in  this  formulation  as  we  are 
currently  developing  a  weak  Lagrange-Galerkin  shal¬ 
low  water  model  on  the  sphere  using  Cartesian  coor¬ 
dinates  [see  Giraldo  (1997)  for  details  of  the  weak  La¬ 
grange-Galerkin  method].  Lor  this  reason,  this  paper 
only  deals  with  the  determination  of  the  departure  points 
in  Cartesian  space. 


Corresponding  author  address:  Francis  X.  Giraldo,  Naval  Research 
Laboratory,  7  Grace  Hopper  Ave.,  Monterey,  CA  93943. 

E-mail:  giraldo @ nrlmry.navy.mil 


McGregor  (1993)  introduced  an  economical  departure 
point  calculation  procedure  that  does  not  involve  inter¬ 
polation.  His  scheme  is  very  efficient  and  accurate  but 
the  implementation  of  the  method  was  only  illustrated 
for  rectangular  grids.  This  paper  shows  the  implemen¬ 
tation  of  McGregor’s  approach  on  unstructured  trian¬ 
gular  grids  such  as  those  composing  the  spherical  geo¬ 
desic  grids. 


2.  Governing  equation 

Let  us  consider  a  simple  equation  for  our  study  such 
that  we  have  analytic  values  for  the  solution  and  the 
trajectories.  In  spherical  coordinates  the  advection  equa¬ 
tion  for  the  variable  ip  is 


U  dip 
a  cos  6  9\ 


v  dtp 
a  30 


+ 


1  9u  1  v 

- 1 - -  tan# 

a  cos 6  9 A  a  99  a 


=  0,  (1) 


where  a  is  the  radius  of  the  sphere,  ( u ,  v)  are  the  zonal 
and  meridional  velocity  components,  and  (A,  6)  are  the 
longitudinal  and  latitudinal  coordinates.  The  first  brack¬ 
eted  term  represents  the  operator  u  •  Vip  and  the  second 
term  represents  ipV  ■  u.  However,  instead  of  using  this 
form,  let  us  look  at  the  Cartesian  form 
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dt 


dip  dip  dip 

U - b  V - b  W - 

dx  dy  dz 


+ 


dll  dv  dw 
^  dx+  dy+  dz 


</>(xA,  t  +  At)ip(xA,  t  +  At)  dPlA 

JaA 


=  0  (2) 

used  in  Giraldo  (1997).  This  equation  can  now  be  writ¬ 
ten  in  the  following  compact  form: 

dip  _ 

V  +  v  •  (up)  =  0,  (3) 

dt 


which  is  the  conservation  form  of  the  advection  equa¬ 
tion  where  ip  is  some  conservation  variable,  and  u  is 
the  velocity  vector.  In  Giraldo  (1997),  this  equation  was 
solved  using  the  weak  Lagrange-Galerkin  method.  Be¬ 
ginning  with  the  method  of  weighted  residuals  with 
weight  (i.e.,  basis)  function  ip 


dip  _ 

- b  V  •  (up) 

dt 


dCl  =  0, 


on  domain  0  and  integrating  by  parts  such  that 

dip  dipip  dip 

ip —  = - ( p — 

dt  dt  dt 

and 


ipV  ■  (up)  =  V  •  (vupip)  —  (up)  •  Vi//, 


we  get 


dipip 

—  +  V  •  (upi//) 
dt 


dip 

pi - b  U  ■ 

1  dt 


Vi p 


dVL  =  0. 


The  basis  functions  ip  are  chosen  such  that  the  second 
term  in  brackets  disappears.  In  other  words,  the  basis 
functions  are  chosen  such  that  they  are  constant  along 
the  characteristics.  This  results  in  the  simplified  system 


dtpip  _ 

- b  V  •  (u< pip)  dD,  =  0 

dt 


along  with 


where 


dip 

—  +  U 

dt 


Vi//  =  0 


dx 

dt 


u(x,  t). 


d_ 

dt 


d 

—  +  u  ■  V 
dt 


(4) 


ip(xD,  t)ip(xD,  t )  dflL 


(5) 


where  A  and  D  denote  arrival  and  departure  points.  [For 
further  details  on  this  method,  refer  to  Giraldo  (1997).] 
However,  this  is  clearly  not  the  only  method  of  solving 
this  equation  using  characteristic-based  methods.  In  the 
case  of  a  divergence-free  flow,  the  Lagrangian  form  of 
(3)  is 


dx 

dt 


u(x,  t). 


(6) 

(7) 


Discretizing  this  equation  by  the  semi-Lagrangian  meth¬ 
od  and  then  applying  the  finite  element  method,  yields 
the  following  relation: 


p(xA,  t  +  At)ip(xA,  t  +  At)  dPlA 

JnA 

=  |  p(xD,  t)ip(xA,  t  +  At)  dClA,  (8) 
JsiA 

which  is  equivalent  to 

p(xA,  t  +  At)  =  ip(xD,  t), 

which  is  the  typical  semi-Lagrangian  formulation  for 
this  equation  regardless  of  which  spatial  discretization 
method  we  select.  Note  that  the  resulting  equations  for 
both  methods  are  equivalent  if  and  only  if  the  flow  is 
divergence  free.  The  main  difference  between  the  two 
approaches  is  that  (5)  depends  on  integration  while  (8) 
on  interpolation.  The  integration  can  be  carried  out  ex¬ 
actly  or  by  Gaussian  quadrature.  The  interpolations,  on 
the  other  hand,  are  a  bit  more  complex  on  unstructured 
grids  and  perhaps  the  best  approach  is  to  use  the  kriging 
method  described  in  Le  Roux  et  al.  (1997),  although 
Lagrange  interpolation  on  the  triangles  could  certainly 
be  used.  However,  the  accuracy  of  both  methods  relies 
mainly  on  how  accurate  the  trajectories  are  calculated; 
namely,  how  best  to  solve  (7).  The  solution  of  this  equa¬ 
tion  on  unstructured  triangular  grids  is  the  scope  of  this 
paper. 


denotes  the  total  (or  Lagrangian)  derivative.  By  virtue 
of  the  Reynolds  transport  theorem,  (4)  now  becomes 

ipip  dPt  =  0, 

which,  after  integrating  along  the  characteristics,  gives 


3.  Test  case 

Numerical  experiments  are  performed  on  the  advec¬ 
tion  equation  on  the  sphere,  which  is  defined  by  (1) 
The  initial  condition  is  given  as  in  Williamson  et  al 
(1992)  by  the  cosine  wave 
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<p(x,  0) 


li 

2 

0 


1  +  cos 


if  r  <  R 
if  r  >  R, 


where  (1  represents  the  element  triangles,  in  a  finite 
element  sense.  In  addition  to  the  L2  norm,  two  more 
measures  are  used,  namely,  the  first  and  second  moments 
of  the  conservation  variable,  which  are  defined  as 


where 
h  =  1, 

r  =  a  arccos[sinf(.  sin0  +  cos0c  cost?  cos(A  —  Ac)], 


ec  =  o. 


2t Ta 

co  =  - 

12  days 


and  the  velocity  field  is  assumed  to  be  constant  and 
given  by 


u  =  +&>(cos0  cosa  +  sin0  cosA  since) 


v  =  —co  sinA  since. 


where  ce  determines  the  axis  of  rotation  of  the  flow  with 
respect  to  the  poles.  As  an  example  a  =  0  yields  flow 
along  the  equator.  The  results  are  reported  for  one  rev¬ 
olution  of  the  initial  wave  that  takes  12  days  to  complete 
one  full  revolution  about  the  sphere.  By  using  the  map¬ 
ping  from  spherical  to  Cartesian  space 


M,  = 


i p(x,  t)  dCl 


<Pexac,(X,  t )  d£l 


and 


M,  = 


i p(x,  t)2  dD, 


<Pexac,(X,  02  d£l 


(12) 


(13) 


These  values  measure  the  conservation  properties  and 
dissipation  of  the  numerical  method,  respectively. 


4.  Trajectory  calculations 


x  =  a  cos  9  cos  A 


The  exact  solution  to  (3)  can  be  written  as 


y  =  a  cos  0  sinA 


<Pex ac,(x,  t  +  At)  =  cp(x  -  Am,  t). 


where 


z  =  a  sind. 


(9) 


A  =  arctanf  — 


6  =  arcsin 


-(:)■ 


(10) 


we  can  write  the  initial  conditions  in  terms  of  Cartesian 
coordinates.  This  results  in  the  following  velocity  field: 

u  =  — u  sinA  —  v  sind  cosA 

v  =  +m  cosA  —  v  sind  sinA 


w  =  +v  cos6. 


along  with  the  following  analytic  solution: 

^exact(X,  t)  =  Cf(X  -  fU,  0), 

which  is  the  solid  body  rotation  of  the  cosine  wave  about 
the  axis  defined  by  a.  Note  that  the  mapping  from  spher¬ 
ical  to  Cartesian  space  is  only  done  once  at  the  begin¬ 
ning  in  order  to  define  the  problem.  From  then  on,  the 
problem  is  solved  in  Cartesian  space.  The  L2  error  norm 
is  defined  as 


[cp(x,  t)  -  <pexact(x,  f)]2  dti 


[<Pexac,(X,  t)]2  dfl 


,  HD 


By  applying  a  rotation  transformation  as  in  McDonald 
and  Bates  (1989),  we  get  the  arrival  points  in  the  rotated 
space 

A^  =  arctan[cos0A  sin(AA  —  AJ, 

cos0A  cos(Aa  -  AJ  cosOa  +  sin#A  sin#J 

0'A  =  arcsin[sin0A  cos0„  —  cos0A  cos(Aa  —  A„)  sin0„], 

which  now  consist  of  motion  about  the  equator,  only. 
Note  that  since  all  the  points  are  stored  in  Cartesian 
space  we  must  use  (10)  in  order  to  get  the  arrival  points 
in  terms  of  spherical  coordinates.  The  departure  points 
in  the  rotated  space  are 

,  ,  Afa> 

a;  =  a;  - 

a 

On  =  01, 

where  Aa  =  0  and  0a  =  a.  Using  the  inverse  transfor¬ 
mation,  we  get 

A»  =  A,,  +  arctan[cos0„  sinAJ,, 

cos cosAq  cos 9a  —  sind^  sindj 

0D  =  arcsin[cos0^,  cosA^  sindn  +  sin0^  cos0a],  (14) 

which  are  the  departure  points  in  the  original  (unrotated) 
spherical  space.  The  departure  points  in  Cartesian  space 
can  now  be  obtained  using  the  mapping  (9).  Equation 
(14)  gives  us  the  exact  trajectories  and  so  we  must  devise 
a  scheme  that  best  approximates  this  solution. 
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a.  Ritchie’s  method 


There  are  many  ways  of  integrating  (7)  but  the  most 
common  form  in  plane  space  has  been  the  midpoint  rule, 
namely 


xM  =  xA 


(15) 


b  = 


a 


a 2  -  A txA  ■  uM  + 


UM 


1/2’ 


which  is  similar  to  the  result  in  Ritchie  (1987)  but  for 
a  two-time-level  scheme.  From  Ritchie  (1987),  we  get 
the  departure  points  from  the  relation 


which  defines  a  recursive  scheme  because  xM  is  given 
implicitly  in  the  relation.  In  addition,  this  scheme  also 
requires  interpolation  because  xM  will  generally  not  fall 
on  a  grid  point.  Usually,  between  three  and  five  iteration 
loops  [see  Ritchie  (1987)]  are  required  to  converge  to 
the  solution  at  which  point,  the  departure  point  is  calcu¬ 
lated  by 


2xi 


-xM  -  xA, 


which  can  be  simplified  by  using  the  definition  of  a  dot 
product 

XM  •  XA  =  IxJ  lxAl  COS0 


to 


X 


D 


-  Aru 


xc  +  xA 
2  cost? 


However,  we  have  said  nothing  about  the  order  of  the 
interpolations  but  obviously  the  higher  the  order  of  the 
interpolant,  the  better  the  trajectory  accuracy,  but  the 
greater  the  computer  time  as  well.  The  midpoint  rule 
yields  second-order  accuracy  and  has  been  used  quite 
successfully  in  2D  planar  space. 

On  the  sphere,  the  midpoint  rule  has  to  be  modified 
such  that  the  new  departure  points  computed  by  (15) 
remain  on  the  surface  of  the  sphere.  In  other  words, 
after  each  iteration  we  must  apply  the  projection 

a 

xo  —  :  ;xD, 

lxDl 

where  a  is  the  radius  of  the  sphere. 

In  fact,  Ritchie’s  method  simplifies  to  the  midpoint 
rule  on  the  surface  of  a  sphere.  In  his  method,  we  start 
with  a  Taylor  series  expansion  about  the  midpoint  ( t  + 
At/2 )  up  to  second  order  to  get 

(  A t\  At  dxf  A t\ 

x((  +  Ar)  =  +  — j  +  y  —  jr  +  — j  +  O(At)2. 

Now,  let  xA  =  x(t  +  At)  and  xM  =  x(t  +  At/2)  be  the 
arrival  and  midpoints,  respectively.  After  rearranging, 
we  get 

At 

xm  —  xA  —  yUM, 

which  is  exactly  equal  to  (15),  where  uM  =  u(x  t  + 
At/2).  However,  we  need  to  add  a  correction  factor  b  in 
order  to  constrain  the  new  iterated  point  to  remain  on 
the  sphere.  Thus,  we  require 


If  we  take  the  midpoint  to  be  the  average  between  the 
arrival  and  departure  points  and  then  project  it  on  to  the 
sphere,  we  arrive  at  the  relation 


By  once  again  enforcing  that  the  point  remains  on  the 
sphere,  we  get 


b  = 


1  +  cos2 9 


and  by  using  the  trigonometric  identity  cos2  0  =  2  cos 2  9 
—  1  we  now  get 


b  = 


1 

cos  O' 


which  recovers  Ritchie’s  result.  Thus,  Ritchie’s  method 
is  exactly  the  midpoint  integration  rule  with  a  modifi¬ 
cation  that  projects  the  iterated  point  onto  the  surface 
of  the  sphere. 


b.  McGregor’s  method 

Another  possibility  for  integrating  (7)  is  the  noniter¬ 
ative  scheme  of  McGregor.  From  McGregor  (1993),  we 
write  a  Taylor  series  expansion  for  the  departure  point 
x(f)  about  the  arrival  point  x(t  +  At)  along  the  char¬ 
acteristics  as 

N  ( —  A  An  rln-v 

X(f)  =  x(f  +  At)  +  2 - ; — +  At) 

„=i  n\  at 


such  that 


(16) 

where 


IxJ  =  a. 


Taking  the  magnitude  of  (16)  and  rearranging  we  get 


and 


+  0(At,  Ax)N+1, 


d  a  „  _ 

—  = - b  u  V  =  u  V 

dt  dt 
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u  =  u^x(f  +  At),  t  +  — 

[see  McGregor  (1993)].  Therefore,  the  only  thing  left 
to  do  is  to  obtain  the  derivatives  in  the  gradient  operator 
V.  If  we  are  using  rectangular  grids  then  we  can  write 
the  derivatives  in  the  standard  centered  finite  difference 
form 


dip 

dx 


(x,  y,  z) 


<p(x  +  Ax,  y,  z)  ~  <p{x  —  Ax,  y,  z) 
2Ax 


(see  the  appendix  for  the  proofs  and  derivations  con¬ 
cerning  these  basis  functions).  This  is  a  necessary  con¬ 
dition  for  a  consistent  and  monotonic  interpolation. 
These  natural  coordinates  can  now  be  used  as  the  finite 
element  basis  functions.  Integration  by  parts  reveals  that 
any  integral  of  the  following  form  involving  these  basis 
functions  can  be  obtained  in  closed  form  by  the  relation 


dil 


deter  Aal/3ly\ 

(a  +  /3  +  y  +  2)! 


(19) 


dip 

dy 


(x,  y,  z) 


dip 

dz 


(x,  y,  z) 


cp(x,  y  +  Ay,  z)  ~  <p(x,  y  -  Ay,  z) 
2  Ay 

<p(x,  y,  z  +  Az)  -  <p(x,  y,  z  -  Az) 
2A  z 


yielding  a  second-order  accurate  derivative  in  space. 
Higher  derivatives  are  also  readily  available  reapplying 
this  relation.  But  what  if  the  grid  is  unstructured?  In 
this  case,  we  must  compute  the  derivatives  in  a  finite 
element  sense.  But  first,  let  us  introduce  the  linear  tri¬ 
angular  finite  element  basis  functions  on  the  sphere. 


5.  Basis  functions 

Linear  natural  coordinates  on  a  triangle  in  3D  Car¬ 
tesian  space  can  be  written  as 


' h(x,  y<  z)  = 


a, x  +  b,y  +  CjZ 

deter  A 


where 


and 


(17) 


a,  =  yjzk  -  ykZj,  bt  =  xkZj  -  Xjzk, 
c,  =  Xjyk  -  xkyp 


This  relation  is  almost  identical  to  the  closed  form  so¬ 
lution  for  linear  triangles  on  the  plane  given  by  Silvester 
(1969).  For  the  special  case  that  the  three-dimensional 
domain  lies  entirely  on  a  plane,  both  integration  rules 
are  equivalent.  Using  these  basis  functions,  we  can  now 
construct  the  derivatives  at  the  grid  points  in  a  finite 
element  sense. 


Derivatives 

We  can  construct  the  derivatives  in  the  following 
manner:  let 


3 


(pe(x,  y,  z)  =  Yj  ^jix,  y,  z)<Pj 

7=1 


(20) 


denote  the  value  of  <pe  within  the  element,  i//,  the  basis 
functions,  and  <pj  the  value  of  the  conservation  variable 
at  the  vertices  (grid  points)  of  the  element  in  question. 
From  (20)  we  get  the  derivatives  within  the  element  to  be 


dip” 

dx 


(x,  y,  z)  = 


1 


difa 

—(x,  y,  z)<Pj 


*i 

x2 

X2 

deter  A  = 

y3 

y3 

Zi 

z2 

where  i,  j,  k  are  cyclical,  that  is,  if  i  =  1,  then  j  =  2, 
and  k  =  3,  and  so  on.  By  using  the  definition  of  the 
natural  coordinates  (17)  and  the  fact  that  the  three  nodes 
on  each  triangle  define  a  plane 

N  •  (x  -  x,)  =  0, 

where  N  is  the  outward  pointing  normal  to  the  triangle 
and  defined  by 

N  =  (x2  —  x,)  X  (x3  -  x,) 


k 

z2  -  z1 

Z3  -  Zi 


(18) 


«  y 2  -  Ti 

-  *i  y3  ~  yi 

it  can  be  shown  that  the  natural  coordinates  satisfy  the 
condition 


fa(x,  y,  z)  +  \jj2{x,  y,  z)  +  ijJiix,  y,  z)  =  1 


d(pe 

dy 


U,  dif/j 

(x,  y,  z)  =  Zj  —(x,  y,  z)tPj 


j=  i 


dy 


d(pe 

dz 


(x,  y,  z)  = 


1 


j=  i 


diffj 

— (x,  y,  z)<pp 
dz 


(21) 


where 


dx 


(x,  y,  z) 


d</t 

—(X,  y,  z) 


dz 


(x,  y,  z) 


deter  A 

b; 

deter  A 


deter  A ' 


from  (17). 

However,  we  need  the  derivatives  on  the  grid  points 
and  not  within  the  elements.  If  we  knew  these  deriva¬ 
tives,  then  we  could  write  them  as 
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dcpe  ^  dip, 

-T~(x,  y,  z)  =  2  y>  z)— 
dx  /=  i  ox 


7=1 

3 


d(pe  -A  Sty 

—  (X,  y,  z)  =  2  lM*.  y.  z)  — 

dy  .7=1  3y 


7=1 

3 


d2(pe  ^  32<p. 

— <*  7.  D  =  2  *M.  y. 


7=1 

3 


d2cpe  ^  d2<p- 

— (x,  y,  z)  =  2  <M*.  y>  ti— 

ay2  y=i  ay2 


9(p‘  ^  dipj 

— (x,  y,  z)  =  2  ’A/*-  y~  z)  ~T-  (22) 

dz  7=1  3z 


7=1 

3 


av  ^  a2^p. 

— (x,  y,  z)  =  2  </';(*.  y>  z)tt- 

dz2  7=1  dz2 


Equating  these  relations  with  (21)  and  employing  the 
finite  element  method,  we  can  construct  a  set  of  integral 
equations,  namely 


dtp,  f  di/t. 
dtl—  =  ip—  dVL<pj 


dx 


dx 


dq>j  f  dl p: 

^ipj  dil—  =  ip—  d£hpj 

Jn  ft 

dtp :  f  dip j 

ip/ipj  di\- —  =  ip , - dfi,(Pj, 


dz 


dz 


(23) 


for  i,  j  =  1,  .  .  .  ,  3  which  is  symmetric  and,  depending 
on  the  node  number  ordering,  tightly  banded.  Using  Eqs. 
(17)  and  (19),  Eq.  (23)  results  in  the  following  elemental 
matrix  equations: 


deter  A 


24 


2  1  1 
1  2  1 
1  1  2 


deter  A 
24  deter  A 


d<Pi 
dx 

dcp  i 

dx 
d<P3 
dx 

a,<p,  +  a2(p2  +  a3<p3 

a,<p,  +  a2cp2  +  a3cp3  • 

«i<Pi  +  a2<P2  +  a3<P3 


(24) 


that  is,  assuming  we  knew  the  second  derivatives  at  the 
grid  points.  Taking  the  derivatives  of  (22)  we  get 

3V.  v  d^f  sd<pJ 

— (x,  y,  z)  =  2j  —(x>  y-  ti— 

OX2  j=  1  ox  ox 

d2(pe  ^  diPj  dcpj 

— (x,  y,  z)  =  2j  —(x>  y-  ti¬ 
dy-  7=1  dy  dy 

aV  2  dip  d<pj 

— (x,  y,  z)  =  2j  ~(x’  y>  z)t-. 

OZ2  7=1  dz  dz 

where  the  derivatives  of  ipj  are  the  first  derivatives  at 
the  grid  points  obtained  from  (24).  Equating  these  two 
sets  of  relations  and  once  again  employing  the  finite 
element  method,  we  obtain  the  equations 

d2(p:  f  dipj  d<p; 

iptf,  dVL-^  =  ip,—1  d£l— 

a  dx  Jn  dx  dx 

d2(p,  f  dlp;  dipj 

^  d[l — ~r  =  ip,—  d£l  — 

1  dy2  Jj/'dy  dy 

d2i Pj  f  dipj  d<Pj 

ipjipj  =  <//,—  dtt-^,  (25) 

dz2  Jn  dz  dz 

which  simplify  to  the  following  element  matrix  relations 

d2<Pi 


with  similar  relations  for  the  y  and  z  terms.  The  global 
system  composed  of  these  elemental  matrix  equations 
requires  the  inversion  of  a  sparse  but  tightly  banded 
matrix.  This  potential  bottleneck  can  be  bypassed  by 
diagonalizing  the  elemental  equations  yielding 


deter  A 


24 


2  1  1 
1  2  1 
1  1  2 


dx2 

d2<p2 


dx2 

d2<?3 

dx2 


1  0  0 
0  1  0 
0  0  1 


d<Pi 

dx 

dy2 

dx 

dy3 

dx 


1 


4  deter  A 


a,<P!  +  a2<p2  +  a3(p3 

ai<Pi  +  a2(p2  +  a3(p} 

a,<Pi  +  a2cp2  +  a3(p3 


deter  A 


24  deter  A 


d(P  i 
1  dx 
diP  i 


+  a 


dip2 
1  dx 
d<P2 


+  a 


dip} 

dx 

d(p3 


a,—  +  a2 — -  +  a, — - 
dx  dx  '  dx 


dip i 

a , - 1-  a? 

dx 


dip 2 

dx 


+  a 


dip3 
’  dx 


(26) 


which  can  be  solved  explicitly  for  the  derivatives.  The 
second  derivatives  within  the  element  can  now  be  writ¬ 
ten  as 


where  once  again,  the  relations  for  the  y  and  z  terms 
are  immediately  obvious.  The  diagonalized  version 
yields 
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8  9  10 

Fig.  1.  The  element-wise  contribution  of  the  surrounding  elements 
to  the  grid  points. 


1  0 
0  1 
0  0 


0 

0 

1 


a2Pi 

dx2 
d2<P 2 
dx 2 
d2<p} 
dx 2 


1 


4  deter  A 


d<Pt  ,  d<p2  dcp3 

a  j - I-  a2 - 1-  a3 — 

dx  dx  dx 

d<Pi  ,  d(Pi  ,  d(p3 

ax - r  a2 - r  a 3 - 

dx  dx  dx 

d<pj  d(p2  dip3 


dx 


+  ci2 — -  +  a 
dx 


dx 


which,  again,  does  not  require  the  inversion  of  a  large 
global  matrix.  Higher-order  derivatives  can  be  obtained 
by  continuing  this  process.  Note  that  Eqs.  (24)  and  (26) 
may  appear  to  imply  that  dtpjdx  =  dcpjdx  =  d(p3/dx 
and  similarly  for  all  derivatives.  However,  these  are  el¬ 
emental  equations  and  the  global  equations  are  obtained 
by  summing  the  contribution  of  all  the  triangular  ele¬ 
ments  surrounding  each  node  point.  Because  the  basis 
functions  are  linear,  then  the  derivatives  within  the  el¬ 
ement  will  be  constants  so  that  Eqs.  (24)  and  (26)  say 
that  the  contribution  of  each  triangular  element  to  its 
three  vertices  is  a  constant  value.  In  Fig.  1  we  can  see 
that  grid  point  5  will  get  contributions  from  elements 
1-6,  while  grid  point  6  will  get  contributions  from  el¬ 
ements  5-10.  As  a  result,  all  of  the  gridpoint  derivatives 
will  have  unique  values. 


6.  Derivative  accuracy 

Before  we  study  the  accuracy  of  the  trajectory  cal¬ 
culations  using  the  results  from  the  previous  two  sec¬ 
tions  it  is  wise  to  first  know  the  accuracy  of  the  deriv¬ 
atives  themselves.  Let  us  apply  the  3D  linear  basis  func¬ 
tions  on  a  plane.  We  perform  these  numerical  experi¬ 
ments  on  the  plane  using  a  structured  grid  in  order  to 
compare  our  finite  element  derivatives  with  the  standard 
centered  finite  differences.  On  this  plane,  let  z  =  c  where 
c  is  some  constant.  Therefore  our  basis  functions  sim¬ 
plify  as  follows: 

=  a-jX  +  b,y  +  cf 

deter  A 


Table  1.  Derivative  accuracy  using  centered  finite  differences  for 
the  cosine  function  in  planar  space  on  the  crisscross  grid. 


Sf 

Sf 

92f 

d2f 

92f 

Grid  points 

dx 

L2 

dy 

L2 

dx2 

f-‘2 

dy2 

L2 

dxdy 

it 

X  11 

0.2498 

0.2498 

0.5076 

0.5076 

0.5174 

21 

X  21 

0.1127 

0.1127 

0.3655 

0.3655 

0.3288 

41 

X  41 

0.0346 

0.0346 

0.2146 

0.2146 

0.2113 

61 

X  61 

0.0174 

0.0174 

0.1602 

0.1602 

0.1667 

81 

X  81 

0.0108 

0.0108 

0.1325 

0.1325 

0.1430 

h  = 

=  yj  ~ 

y*. 

ja¬ 

il 

1 

Xp  c, 

=  Xjyk 

-  xkyp 

and 


deter  A  = 


*1 

1 


x2  x3 

y2  y3 

l  l 


where  z  =  c  has  been  factored  out  from  all  of  the  re¬ 
lations.  Recall  that  these  functions  have  the  exact  in¬ 
tegration  rules  obtained  by 


Ip"  1^2  dil 


deter  Aalfilyl 
(a  +  (3  +  y  +  2)!  ’ 


which  now  yields  the  usual  linear  triangular  finite  ele¬ 
ment  basis  functions  on  the  plane  as  in  Silvester  (1969). 
To  test  the  accuracy  of  the  numerical  derivatives,  we 
shall  use  the  cosine  function 


m 


h 

2 

0 


1  +  cos 


if  r  <  R 
if  r  >  R 


with 

r  =  V(x  —  xc)2  +  (y  -  yc)2. 


(xc,  yc)  =  (0,  0), 


[x,  y]  e  [-1,  +1],  R  = 


By  defining  a  normalized  L2  norm  such  as 


df 

1 

Ja 

,  (x)  ,  (x) 

dx  dxexM 

2 

dtl 

dx 

f 

Ja 

df  f  t 
to  (X) 

L,*A'exact 

2 

dn 

for  all  of  the  derivatives,  we  can  now  measure  the  ac¬ 
curacy  of  our  numerical  derivatives.  Tables  1,  2,  and  3 
list  the  first  and  second  derivative  results  for  various 
grid  sizes  using  centered  finite  differences,  finite  ele¬ 
ments  with  a  full  matrix,  and  finite  elements  with  a 
diagonalized  matrix,  respectively.  An  example  of  the 
structure  of  the  grids  used  is  illustrated  in  Fig.  2  for  the 
11X11  point  case.  The  tabulated  results  show  that  the 
finite  element  derivatives  are  superior  to  the  finite  dif¬ 
ferences  but  this  approach  requires  the  inversion  of  a 


where 
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Table  2.  Derivative  accuracy  using  finite  elements  with  a  full 
matrix  for  the  cosine  function  in  planar  space  on  the  crisscross  grid. 


Sf 

3f 

s2f 

32f 

d2f 

Grid  points 

dx 

L2 

dy 

L2 

dx 2 

L2 

fy 2 

^-2 

dxdy 

ii 

X 

11 

0.1536 

0.1536 

0.5678 

0.5678 

0.3170 

21 

X 

21 

0.1047 

0.1047 

0.3806 

0.3806 

0.2959 

41 

X  41 

0.0262 

0.0262 

0.2251 

0.2251 

0.1886 

61 

X 

61 

0.0162 

0.0162 

0.1712 

0.1712 

0.1508 

81 

X  81 

0.0087 

0.0087 

0.1422 

0.1422 

0.1263 

matrix.  For  this  reason,  we  have  also  included  the  di¬ 
agonalized  version,  which  does  not  require  this  inver¬ 
sion.  The  results  show  that  the  diagonalized  version  is 
inferior  to  the  full  matrix  version  but  if  we  are  concerned 
with  efficiency,  then  this  faster  version  is  more  appro¬ 
priate.  Let  us  now  look  at  the  formal  order  of  accuracy 
analysis  of  the  finite  element  type  numerical  derivatives. 


Order  of  accuracy  analysis 

The  structured  crisscross  grid  illustrated  in  Fig.  2  was 
selected  for  our  study  because  it  has  the  least  amount 
of  biasing.  This  is  important  because  the  spherical  geo¬ 
desic  grid  has  little  or  no  biasing  on  the  sphere.  The 
crisscross  grid  has  two  types  of  gridpoint  contributions 
from  the  surrounding  elements  that  are  illustrated  in 
Figs.  3  and  4.  The  first  type  of  contribution  (illustrated 
in  Fig.  3)  has  derivatives  given  by 


dfu 

dx 


C/h 


2Ax 


<9(Ax2) 


Vfu 

dx2 


(fi+2. 


~  2fi.j  + 

4Ax2 


+  0(Ax2) 


Fig.  2.  The  11X11  point  crisscross  grid  used  for  the  derivative 
analysis. 


tered  finite  differencing.  The  second  type  of  contribution 
(Fig.  4)  has  the  derivatives 

dfij  =  ~  /,- u+i)  +  2 -  y;_w) 

dx  8Ax  8Ax 

+  - lizldzA  +  0(Ax2  Ay2), 

8  Ax  ^ 


Wjj 

dxdy 


(f, 


+  1J+1 


/, 


1J+1 


/h 


1J-1 


+  fi-u- 1) 


4AxAy 


+  0( Ax2,  Ay2), 


where  the  derivatives  in  y  are  immediately  obvious.  The 
orders  O  given  in  these  relations  denote  the  order  of 
accuracy  of  the  derivatives.  These  derivatives  are  all 
second-order  accurate  in  both  x  and  y.  In  fact,  they  yield 
the  exact  same  derivative  formulas  obtained  with  cen- 


Table  3.  Derivative  accuracy  using  finite  elements  with  a  diago¬ 
nalized  matrix  for  the  cosine  function  in  planar  space  on  the  crisscross 
grid. 


3f 

if 

d2f 

32f 

32f 

Grid  points 

dx 

L2 

dy 

L2 

dx2 

L2 

3y 2 

i-2 

dxdy 

11 

X 

11 

0.2879 

0.2879 

0.7452 

0.7452 

0.6645 

21 

X 

21 

0.1341 

0.1341 

0.4519 

0.4519 

0.3874 

41 

X  41 

0.0382 

0.0382 

0.2718 

0.2718 

0.2402 

61 

X 

61 

0.0207 

0.0207 

0.2095 

0.2095 

0.1884 

81 

X  81 

0.0121 

0.0121 

0.1764 

0.1764 

0.1618 

i,  j+1 


i+1,  j 


Fig.  3.  The  first  type  of  element  contribution  to  the  grid  points  of 
the  crisscross  grid. 
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i-1 ,  j-1  i,  j-1  i+1,j-1 


Fig.  4.  The  second  type  of  element  contribution  to  the  grid  points 
of  the  crisscross  grid. 


d2f,j  =  c fi+2,j+2  -  2 fu+2  +  /;_2J+2) 
dx 2  64Ax2 

4 (fi+2J+l  ~  2fiJ+l  +  /,_2J+1) 

64Ax2 

6(f,+2j  -  2f,j  +  /,_2,,) 

64Ax2 

Ww-i  -  2/,,,-!  +  ,/; ,  ,  ,i 

64Ax2 

,  (fi+Xj-2  ~  Vi.j-2  +  fi-2,j-2)  ,  ,  A  ,s 

+ - ^ - +  OAr,  Ar). 

and 


^~/ij  _  (fi+2,j+2  fi+2,j-2  fi-2,j+2  4  fi-ZJ-f) 

dxdy  64AxAy 

2(/j+2J+l  ~~  Ji+2J—1  ~~  fi-2J+ 1  4"  J i—2,j—  l) 

64AxAy 

2(/j+lj+2  _  fi+lJ—2  —  fi  !,_/+ 2  5"  ft-lJ-2 ) 

64AxAy 

4(/i+ij+i  —  ,/,-h  j-i  —  f i  i  ,y+ 1  4“  ft- Ur  i) 

64AxAy 

+  0(Ax2,  Ay2), 

which  are  all  second  order.  Thus  even  with  the  diago¬ 
nalized  version  of  the  derivatives,  we  are  guaranteed  an 
order  of  accuracy  similar  to  centered  finite  differences, 
regardless  of  the  structure  of  the  grid.  This  order  of 
accuracy  analysis  is  given  only  to  compare  the  finite 
element  derivatives  to  the  finite  difference  derivatives. 


Fig.  5.  The  differencing  stencil  for  the  finite  element  derivatives 
on  the  spherical  geodesic  grid. 


However,  the  real  interest  is  in  determining  the  finite 
element  derivatives  on  the  hexagonal-type  stencils  that 
arise  in  the  spherical  geodesic  grids  and  is  illustrated  in 
Fig.  5.  In  this  case,  we  arrive  at  the  following  deriva¬ 
tives: 

d/ij  _  */;+(i/2)j+i  ~  fi-(i/2\j+i  )  +  2(/,+1,  - 

dx  6Ax  6Ax 

,  (/i+(i/2)j-i  —  /i-(l/2)j-i)  ,  ,  A  ^ 

+  - — -  +  <9(Ax2,  Ay-), 

6Ax 

4//J  (/;+(l/2W+l  ~  fi  I  ( 1/2)./  l)  (fi-(U2),j+ 1  —  f i-(U2),j-\) 


dy  4Ay 

+  <9(Ax2,  Ay2), 

2  -  2/;,.,  +  ./;  w,,) 
dx2 


4Ay 


+ 


+ 


+ 


+ 


36Ax2 

4(/i+(3/2)j+l  _  /i+(l/2)J+l  1/2), 7+ 1  4/',  (3/2)./ I  l) 

36Ax2 

2(2 f+2J  +  /f+ij  -  6/,y  +  +  2/;_,J) 

36Ax2 

4(/;+(3/2)j-1  _ /i+(l/2),j-l  f i-(\/2),j- 1  4/)  (3/2) j  i) 

36Ax2 

(fi+lJ-2  ~  lfi.j-2  +  f i— 1,7—2) 


36Ax2 

a2/,.,  (/,+lJ+2  -  2/,+1J  +  /,+1J_2) 

dy2  16  Ay 2 

2(/,.2  -  2/y  + 


+  0(Ax2,  Ay2), 


+ 


+ 


16  Ay 2 

(f-ij+2  ~  2/;._1J  +  fi-hj-2) 

16  Ay 2 


+  <9(Ax2,  Ay2), 
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Table  4.  Trajectory  accuracy  for  the  exact  trajectories,  midpoint 
rule,  and  McGregor’s  noninterpolating  scheme  for  values  of  N  =  1 , 
.  .  .  ,  4  with  a  =  1 . 1 3  on  the  spherical  geodesic  grid  with  642  points 
after  one  revolution  for  a  =  0. 


Trajectory  method 

IWL 

IML 

Mi 

M, 

Exact 

0.0000 

0.0917 

1.0071 

0.9783 

Midpoint  rule 

0.0026 

0.1132 

1.0069 

9.9809 

N  =  1 

0.0358 

0.4743 

1.0127 

0.6064 

N  =2 

0.0041 

0.1724 

1.0064 

0.9764 

N  =  3 

0.0003 

0.0918 

1.0069 

0.9854 

N  =  4 

0.0002 

0.0917 

1.0069 

0.9844 

Table  5.  Trajectory  accuracy  for  the  exact  trajectories,  midpoint 
rule,  and  McGregor’s  noninterpolating  scheme  for  values  of  N  =  1 , 
.  .  .  ,  4  with  a  =  2.27  on  the  spherical  geodesic  grid  with  2562  points 
after  one  revolution  for  a  =  0. 


Trajectory  method 

IWL 

IML 

M, 

m2 

Exact 

0.0000 

0.0195 

0.9988 

0.9996 

Midpoint  rule 

0.0008 

0.0386 

0.9988 

0.9995 

N  =  1 

0.0358 

0.4684 

1.0050 

0.6271 

N  =  2 

0.0041 

0.1506 

0.9980 

1.0001 

N  =  3 

0.0003 

0.0210 

0.9987 

1.0057 

N  =  4 

0.0002 

0.0206 

0.9987 

1.0057 

and 


d  jij  _  (fi+lji-2  fi  '  I ,j  2  Ji-l,j+2  Ji-lJ-l) 

dxdy 


+ 


+ 


24AjrAy 

2(/i+(3/2)j+l  —  /;+(3/2)J-l  .//— C3/2).y+l  /i'-(3/2),j-l) 

24AxAy 

2(/i+(l/2X/+l  —  fi-V(V2),j  I  ~  fi  (I/2),J!|  fi-i\n)J  l) 


24AaA  y 

+  0(Ax2,  Ay2), 

which  are  also  second-order  accurate. 


7.  Trajectory  accuracy 


In  a  similar  manner  described  in  McGregor  (1993), 
we  define  the  accuracy  of  the  trajectories  by  the  nor¬ 
malized  L2  norm 


[xD  -  x“ac']2  d£l 


[xD  -  xA]2  dn 


(27) 


The  results  for  various  methods  of  obtaining  the  trajec¬ 
tories  are  illustrated  in  Tables  4  and  5  where  the  spher¬ 
ical  geodesic  grid  contains  642  points  with  a  Courant 
number  cr  =  1.13,  and  2562  points  with  a  =  2.27, 
respectively.  These  results  are  all  shown  for  a  =  0 
meaning  that  the  wave  moves  along  the  equator.  A  sche¬ 
matic  of  the  grid  containing  2562  grid  points  is  illus¬ 
trated  in  Fig.  6.  The  results  point  toward  the  same  con¬ 
clusions,  namely,  that  McGregor’s  scheme  is  extremely 
good  and  that  it  increases  in  accuracy  as  the  number  of 
terms  in  the  Taylor  series  N  is  increased.  However,  very 
little  is  gained  beyond  values  of  4,  which  is  in  agreement 
with  the  findings  in  McGregor  (1993).  For  this  reason, 
results  for  N  >  4  are  not  shown. 


8.  Conclusions 

The  determination  of  departure  points  are  explored 
for  spherical  geodesic  grids  in  Cartesian  space.  The  mid¬ 


point  rule,  which  is  an  interpolating  and  iterative 
scheme,  is  compared  against  McGregor’s  noninterpo¬ 
lating  and  noniterative  method.  McGregor’s  method 
yields  better  results  but  no  benefits  are  gained  by  using 
more  than  four  terms  in  the  Lagrangian  Taylor  series 
expansion.  McGregor  (1993)  showed  how  to  apply  this 
scheme  on  rectangular  grids.  This  paper  extends 
McGregor’s  method  to  unstructured  triangular  grids. 
The  difficulty  in  applying  McGregor’s  method  to  un¬ 
structured  grids  is  that  derivatives  at  the  grid  points  need 
to  be  obtained  in  order  to  get  the  higher-order  Taylor 
series  expansion  terms.  This  can  be  done  for  unstruc¬ 
tured  grids  by  constructing  the  derivatives  in  an  element 
by  element  approach.  This  approach  is  illustrated  for 
the  unstructured  triangular  grids  composing  the  spher¬ 
ical  geodesic  grids  by  using  the  linear  triangular  basis 
functions  introduced  in  Giraldo  (1997).  Once  these 
functions  have  been  defined,  we  can  then  apply  the  strat¬ 
egy  for  forming  derivatives  illustrated  here.  The  nu¬ 
merical  derivatives  are  compared  against  analytic  so¬ 
lutions  for  the  cosine  hill  function  and  its  derivatives. 


Fig.  6.  The  spherical  geodesic  grid  with  2562  points. 
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These  results  show  that  the  numerical  derivatives  are 
quite  accurate  especially  for  the  low-order  derivatives. 
An  order  of  accuracy  analysis  is  performed  that  dem¬ 
onstrates  the  order  of  accuracy  of  this  strategy  to  be 
second  order.  Therefore,  it  is  quite  similar  to  the  cen¬ 
tered  finite  difference  approach  used  in  McGregor 
(1993),  but  for  unstructured  triangular  grids. 
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APPENDIX 

3D  Triangular  Basis  Function 


iAi(x) 


(x) 


X  x2  x3 

Xj  X  x3 

y  y3 

u  y  y3 

z  z2  z3 

Zi  z  z3 

lp2(x)  = 


*1 

x2 

x3 

X, 

x2 

x3 

u 

y2 

yi 

y2 

J3 

Zi 

Z2 

Z3 

Zi 

Z2 

23 

Xi 

X2 

X 

yi 

y2 

y 

z  1 

z2 

z 

X, 

X2 

x3 

y  i 

y2 

y3 

Z  i 

z2 

Z3 

Theorem  1.  The  basis  functions 

dfX  +  b:y  +  CjZ 


<Mx)  = 


deter  A 


(Al) 


where 


«.■  =  yjzk  -  ykZj ,  b,  =  Xj 

'j~ 

XjZk, 

and 

x, 

x2 

x3 

deter  A  = 

y3 

y3 

Zi 

Z2 

Z3 

which  can  be  written  more  compactly  in  terms  of  scalar 
and  vector  products  as 


«Ai(x) 


x  •  (x2  X  x3) 
x,  •  (x2  X  x3)’ 


* A2(x) 


x  •  (xt  X  x3) 
x,  ■  (x2  X  x3)’ 


I A3(x) 


X'(x,  X  x2) 
x,  •  (x2  X  x3)’ 


and  finally  we  can  write  (A3)  as 


<Mx) 


x  •  (x2  X  xt) 
x,  ■  (x,.  X  x*)’ 


where  the  following  identities  have  been  used: 


(A3) 


(A4) 


define  a  set  of  linear  cardinal  basis  functions,  that  is. 


<Mx/) 


1 1  if  i  =  j 
[0  if  i  j. 


(x,  X  Xj)  =  -(XjX  x;)  and 
x,  •  (x;  X  xk)  =  Xj  •  ( xk  X  x,.)  =  xk  •  (x,.  X  x;). 
From  (A31) 


Proof.  The  conditions  to  be  satisfied  by  a  linear  in¬ 
terpolation  on  a  triangle  in  three-dimensional  space  are 
the  following: 

X  =  <pfx)xl  +  l j/2(x)x2  +  f//,(x)x3 

y  =  *Ai(x)y,  +  iA2(x)y2  +  <A3(x)y3 

z  =  <Mx)m  +  ' A2(x)z2  +  </'3(x)z3,  (A2) 


<A,(x,) 


x,  •  (x,  X  xk) 
x,  •  (x,  X  xj 


<A,(x2) 


Xj  •  (x,  X  Xt)  _  X<;  •  (Xj.  X  Xj) 
X,  '  (X;  X  Xj)  X,-  •  (Xj.  X  Xj) 


<Ai(x*) 


Xt  •  (Xj  x  Xj)  _  Xj  •  (Xj  X  Xj) 
X;  ■  (Xj.  X  Xj)  Xj  •  (Xj  X  Xj) 


which  just  says  that  the  coordinates  within  the  triangular 
element  are  dependent  on  the  vertices  of  that  element. 
Equation  (A2)  can  be  written  in  the  following  matrix 
form: 


Xj  x2  x3 

iAi(x) 

X 

Ti  y2  y3 

i/k(x) 

= 

y 

_Zj  z2  z3_ 

L  <A3(x)_ 

_z_ 

Using  Cramer’s  rule  to  invert  this  matrix  system  yields 
the  following  relations  for  the  natural  coordinates. 


Q.E.D. 

Theorem  2.  Furthermore,  the  basis  functions  (Al)  sat¬ 
isfy  the  relation  for  a  monotonic  interpolant 

3 

2  ^(x)  =  1  Vx  e  T123, 

i=1 

where  Tj  2  3  represents  the  triangle  composed  of  the  ver¬ 
tices  (x,,  x2,  x3). 

Proof.  Taking  the  sum  of  the  basis  functions  in  (A3) 
and  using  the  identity  (x,  X  x,)  =  0  gives 
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2  <A/(x) 

1=1 

=  X-  [(x,  X  x3)  -  (x2  X  X,)  -  (Xi  X  x3)  +  (Xi  X  X,)] 

x,-(x2Xx3) 

(A5) 

The  numerator  of  (A5)  can  now  be  factored  to  yield 

V  ,  ,  ,  x  [(x2  -  x,)  X  (x3  -  x,)] 

j=i  x3  •  (x2  X  x3) 

The  denominator  can  be  handled  in  a  similar  fashion  to 
yield 

2  x  •  [(x,  —  x,)  X  (x,  —  x,)l 

Y  i k(x)  =  - — - - - — - — ,  (A6) 

£  V  xt  •  [(x2  -  Xl)  X  (x3  -  Xl)] 

where  the  identity 

x;  •  (Xj  X  X,)  =  Xj  ■  (Xj  X  x,-)  =  0 

has  been  used.  The  terms  inside  the  brackets  of  (A6) 
are  exactly  the  components  of  the  normal  vector  to  the 
triangle  (x,,  x2,  x3).  In  other  words 

N  =  (x2  —  x,)  X  (x3  -  x,). 

From  the  definition  of  the  plane  defined  by  this  triangle 
which  is 

N  ■  (x  —  Xj)  =  0 

we  get 

x  •  N  =  x,  •  N, 
which  gives  for  (A6) 

2  x  •  N 

§  =  =  *' 


l//3  =  1  -  i//,  -  lp2 

and  so  r/r3  =1  —  ^  —  7?.  Therefore,  the  integral  becomes 


i//"i//2  i//J  d£l 


Cl(x,y,z) 


=  deter  A 


i-e 


77^(1  -  £  -  P)y  dr] 


d£. 


which  gives  for  the  first  integral  in  brackets  by  virtue 
of  integration  by  parts 


i-f 


17^(1  -  £  -  rf)y  dr] 


V- 


(P  +  1)03  +  2)  •  •  •  (/3  +  y  +  1) 


(1  -  g)o+n-i, 


Q.E.D. 


but  by  completing  the  factorial  in  the  denominate^  we 
get 

f  Rlyl 

77^(1  -  £  -  r])r  dr]  = - — - (1  -  £)f>+r+  K 

,0  (P  +  y+  l )! 

Integrating  the  remaing  terms,  we  get 
i//1“i/>2  i//3  dO, 

deter  Af3lyl(f3  +  y  +  1)! 

(f3  +  y  +  l)!(a  +  l)(a  +  2)  •  •  •  (a  +  j3  +  y  +  2)’ 

and  simplifying  and  completing  the  factorial  yields 

„  ^  deter  Aa!G!y! 

dTl  =  - . 

n(w)  (a+  (3+  y  +  2)\ 

Q.E.D. 


nuy.z) 


Theorem  3.  The  integral  of  any  combination  of  the 
basis  functions  i/r(x)  within  each  triangle  can  be  given 
in  closed  form  as  follows: 


<P?<p2  dil  = 


deter  Aa!/3!y! 

(a  +  (3  +  y  +  2)!  ’ 


(A7) 


where 


v, 

x2 

*3 

deter  A  = 

3b 

y 2 

3b 

M 

Z2 

z3 

Proof.  From  (A2),  the  mapping  from  (x,  y,  z) 
i/j2,  i/j 3)  is 


(A8) 

((/»!, 


dx 

X, 

x2 

x3 

dij/l 

dy 

= 

3b 

3b 

3b 

diji2 

dz 

z, 

^2 

^3_ 

dif 

Let  £  =  i//j,  17  =  i//2  and  from  theorem  2, 
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